Randomly Distributed Delayed Communication and Coherent Swarm 

Patterns 



Brandon Lindley and Luis Mier-y-Teran-Romero and Ira B. Schwartz 



Abstract — Previously we showed how delay communication 
between globally coupled self-propelled agents causes new 
spatio-temporal patterns to arise when the delay coupling is 
fixed among all agents [1]. In this paper, we show how discrete, 
randomly distributed delays affect the dynamical patterns. In 
particular, we investigate how the standard deviation of the time 
delay distribution affects the stability of the different patterns 
as well as the switching probability between coherent states. 

I. INTRODUCTION 

Numerous recent investigations have been devoted to the 
study of interacting multi-agent or swarming systems in 
various natural and engineering fields of study. Investiga- 
tions of interacting systems have revealed the emergence 
of highly complex dynamic behaviors in space and time 
which arise even though the dynamics of a single agent is 
quite simple [2]. In particular, these multi-agent swarms can 
self-organize in complicated spatio-temporal patterns that 
depend on the details of the inter-agent interactions. These 
investigations have been motivated by and had an impact on 
many diverse biological systems such as bacterial colonies, 
schooling fish, flocking birds, swarming locusts, ants, and 
pedestrians [3], [4], [5], [6], [7]. In this paper, we are 
interested in the application that biological analogies have on 
the design of systems of autonomous, inter-communicating 
robotic systems [8], [9], [10], [11] and mobile sensor net- 
works [12]. 

There is great interest to design agent-interaction protocols 
to carry out robotic motion planning, consensus and cooper- 
ative control, and spatio-temporal formation. One methodol- 
ogy is to combine inter-agent potentials with external ones in 
order to achieve multi-agent cooperative motion in a manner 
that is not too sensitive with respect the number of agents. 
Some important applications making use of scalable numbers 
of agents are: obstacle avoidance [10], boundary tracking 
[13], [14], environmental sensing [12], [15], decentralized 
target tracking [16], environmental consensus estimation 
[12], [17] and task allocation [18]. 

Authors have employed very diverse approaches in the 
study of multi-agent systems. Some authors have described 
the swarms at the individual level, writing their models in 
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terms of ordinary differential equations (ODEs) or delay 
differential equations (DDEs) to describe their najectories 
[19], [20], [9]. The addition of noise on the swarm's dynam- 
ics introduces even richer behavior, such as noise-induced 
transitions between different coherent patterns [2], [1]. The 
study of noisy swarm dynamics has benefited from tools from 
statistical physics applied to both first and second order phase 
transitions that have been found in the formation of coherent 
states [21]. 

One important aspect of the understanding and design of 
space-time behavior in communicating robotic systems is 
that of time delay. Time delay arises in latent communication 
between agents, as well as actuation lag times due to inertia. 
Time delays can have interesting and surprising dynamical 
consequences in a system, such as large-scale synchro- 
nization [22], [23], [24], and have been used successfully 
for control purposes [25], [26]. Many of the initial time- 
delay studies focused on the case of one or a few discrete 
time delays. Recently, more complex situations have been 
considered such as the case of having several [27] and 
random time-delays [28], [29]. Another interesting case is 
that of distributed time delays, i.e. when the dynamics of the 
system depends on a continuous interval in its past instead 
of on a discrete instant [30]. 

In the case of swarming systems in stochastic environ- 
ments, it has been observed that the introduction of a discrete 
communication time delay induces a nansition from one 
spatio-temporal pattern to another as the time delay passes a 
certain threshold [1]. It was shown in [1] how the complex 
interplay exists between the attractive coupling and the time 
delay in the transitions between different spatio-temporal 
patterns [31], [32]. Time delays in robotic systems have been 
also studied in the contexts of consensus estimation [17] and 
task allocation [18]; in the latter, the time delays originate 
from the period of time required to switch between different 
tasks. 

In this paper, we consider a swarming model with discrete, 
randomly distributed time delays. We explicitly show how a 
disnibution of delay times perturb the dynamics from the 
single discrete case delay case analytically. We illustrate the 
dynamical effects of delay distributions with varying width 
and show that the system is bistable, and very sensitive to 
choice of initial starting conditions. 

II. Swarm Model 

We investigate the dynamics of a two-dimensional system 
of N identical self-propelling agents that are attracted to each 
other in a symmetric manner. We consider the attraction 



between agents to occur in a time-delayed fashion, due to 
the finite communication speeds and information-processing 
times. Specifically, we focus on the situation in which the 
time-delay is nonuniform across agents: there is one time 
delay for every pair of agents Ty(= T ;! ), for particles ; 
and j. The time delays Ty's are time-independent and are 
drawn independently from a random distribution p T (x). The 
swarm dynamics are described by the following governing 
equations: 



r, =V; 
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v, = (1 - |v,-| 2 ) v,- I £ (r,(r) - rj(t Ty)), (lb) 
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for ; = 1,2..., AT. The position and velocity of the zth agent 
at time t are denoted by r, and v,, respectively. Each 
agent has self-propulsion and frictional drag forces given by 
the expression term (l — |v,| 2 ) v,-. The coupling constant a 
measures the strength of the attraction between agents and 
the communication time delay between particles i and j is 
given by Ty. Note that in the absence of coupling agents tend 
to move in a straight line with unit speed as time tends to 
infinity. 



III. Mean Field Approximation 

We carry out a mean field approximation of the swarming 
system by switching to particle coordinates relative to the 
center of mass and disregarding the noise terms. The center 
of mass of the swarming system is given by 
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We can decompose the position of each particle into 

r,=R + Sr„ * = 1,2...,tV, (3) 
where we'll have 
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Inserting Eq. into the second order system equivalent to 
Eq. (HJ and simplifying we get 
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Summing Eq. over ; and using Eq. (g), we get 
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We now make some approximations on the terms with the 
double sums. For the displacements from the center of mass, 
we have 
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^W^-Jo E 5r 7'( r - T)p T (r)dT = 0, 



(7) 



since Y!j=i ^ r j( f — t) = by Eq. ©. In passing from 
the discrete to the continuous averaging above, we argue 
as follows. The expression YJi=u^j ^ r 7'( f — %) i s me 
average of 8rj(t) at the TV — 1 times t — Ty. Since N 3> 1 
and the times Ty are distributed with density Pt(t), this is 
approximately equal to JJ° 5r/(/ — %)p x {x)dx. 
Similarly, 
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In a purely heuristic manner, we neglect all fluctuation 
terms 5iy(f ) in the dynamics of the center of mass and obtain 
the mean field approximation: 

R = (l-|R| 2 )R- fl ^R(f)-^°°R(r-T)p T (T) £ /T^ . (9) 

where we approximated ^fe! m 1, since we are considering 
large numbers of agents. 

IV. Bifurcations in the Mean Field Equation 

The behavior of the system in the mean field approxima- 
tion in different regions of parameter space may be better 
understood by using bifurcation analysis. This mathematical 
technique will allow us to show how the parameter plane of 
coupling constant a and mean time delay /7 T is divided into 
regions with different dynamical behaviors. 

First we show that Eq. (0 has a uniformly translating 
solution R(f ) = Ro + Vo • t, where Ro and Vo are constant, 
two-dimensional vectors. Inserting the uniformly translating 
state into Eq. (O, we get 

= (1 - |V | 2 ) Vo - aj~ T Pz (T)dT Vo, (10) 



since f{f p T (t)dT = 1. Hence, the speed |Vo| of the uniformly 
translating state must satisfy 
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Tp T (x)dT = 1 -fljU T 
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where ji T is the mean of the p T distribution. We note that 
the direction of motion and starting point Ro are arbitrary. 

The other state of interest is the stationary state R(f ) = Ro, 
for an arbitrary constant vector Ro. In the two-parameter 
space (a,jJ. T ), the hyperbola ajj. T = 1 is in fact a pitchfork 
bifurcation line on which the uniformly translating states are 
born from the stationary state. 

The linear stability of the stationary state is determined by 
the solutions to the characteristic equation of Eq. (0: 



9{X)=a[l- 



-A + A 2 = 0, (12) 



and so involves the Laplace transform of the distribution p T . 

In our numerical simulations of system (|T), we considered 
a truncated Gaussian distribution: 







if t>0 
if T < 0, 
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where jY is the normalization constant. Note that because 
of the truncation, To and Ti are only approximately equal 
to the mean and standard deviation of p T and JV is only 

approximately 1 / ^2nx 2 . 

We approximate the Laplace transform of the truncated 
Gaussian distribution by extending the i ntegra tion range to 
the whole real line and taking jV sa 2m 2 . In addition, 
we approximate the mean and standard deviation of p T as 
jl T w To and (7 T « Ti , respectively. The result is 



p x {x)e- Xx dv, 



(14) 



We use the above approximation to the Laplace transform 
of p T to search for Hopf bifurcation curves in the (a, fl T ) 
plane, by taking A = ico in the characteristic equation (fT2l . 
The equation @(i(D) = is equivalent to: 



a — CO — iC0 = ae 



(15) 



from which we obtain the Hopf curves parameterized by co: 
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In the above expression for li t h(co), the branch of tan in 
(0,7r) should be used, since the complex number on the 
left hand side of Eq. (l5[ is always on the top half plane. 
This family of Hopf curves labeled by re, together with the 
pitchfork bifurcation curve ajj, T = 1 are shown in Figure [TJ 
for various values of ov. 
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Fig. 1. Hopf (blue) and pitchfork (red) branches in a and jl r space. 
The standard deviations of the time-delay distribution p T for the panels 
(a) through (d) are 0, 0.2, 0.4 and 0.6, respectively. Note the change of 
scale in the abscissae. 



When a T = 0, the system exhibits a degenerate point at 
a = 1/2, li z = 2 (Fig. |l(a)[ ), where the Hopf bifurcation 
frequency becomes zero. This is similar, but not equivalent 
to a Bogdanov-Takens bifurcation as is known from previous 
work [31], [32]. Since the point on the Hopf curve in a two- 
parameter bifurcation plane occurs when the Hopf frequency 
becomes zero, we define this point as a Zero Frequency Hopf 
(ZFH) point. 

For (7 T > 0, this ZFH point shifts and a second ZFH point 
appears at a — > °° and /i T — > 0. The location of the two ZFH 
points in the (a, jj, T ) plane is given by {o"zf H , 1 / a zf//)' wnere 



When o~i = 0, the behavior of the mean field in the vicinity 
of the ZFH point is relatively well understood [32], [31] 
and is as follows (see Fig. |l(a)) . In the region between 
the pitchfork and the first member of the Hopf family, the 
stationary state is stable. A simulation of the full system (HJ 
with parameters in this area reveals that indeed the center 
of mass of the agents comes to rest as time progresses 
and the particles spread themselves along a ring with radius 
1 /y/a. Roughly half of the particles move clockwise and the 
other half counterclockwise. Along the first Hopf curve, a 
stable limit cycle is born and the center of mass begins to 
oscillate periodically on a circular orbit. Below the pitchfork 
bifurcation curve aji T = 1, the translating state is stable. 
Finally, we mention that there is a region of bistability in the 
parameter region above the ZFH point (1/2, 2) between the 
pitchfork curve aji z = 1 and the curve aji 2 = 2 (not shown), 
where the center of mass can either translate or rotate. On 
the curve ail 2 = 2 there is a global bifurcation where the 
radius of the orbit diverges and the limit cycle disappears. 

The above discussion helps us understand the bifurcation 
planes in Figs. |l(b)| through |l(d)| Most significantly, we see 
that the parameter region where the stationary state is stable 
decreases in size as the width of the time delay distribution 
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Fig. 2. Two stable attractors for the swarm dynamics. Here a = 2, \l x = 2, 
and o~ T = 0. 15. The number of particles is set to be N = 150. The final state 
shown for both simulations is t = 300. Panel (a) depicts the rotating state 
at three snapshots at times t =297.6, 298.8, 300, in red, green and blue, 
respectively. Panel (b) depicts the ring state. 

widens. Hence the system has a higher tendency to behave 
in an oscillatory manner for wider time delay distributions. 
This effect has been corroborated in numerical simulations 
(results not shown). 

V. Numerical Simulations 

We analyze the dynamics of system ([TJ by solving the sys- 
tem of DDEs numerically. We use Heun's method together 
with quadratic Lagrange interpolation to evaluate the time- 
delayed terms of Eqs. ((TJ. Overall, the numerical method 
is second order with respect to the step-size At. For all 
simulations we take the agents to be uniformly distributed 
in a random fashion within the unit box < x < 1 and 
< y < 1, and each particle is initially at rest v 7 = 0. 
Moreover, since we are interested in investigating the time- 
asymptotic behavior, for all numerical experiments the time 
of integration is long enough to allow transients to decay. 

In [31], [32] it was shown that for the parameter set a = 2, 
T = 2 (fixed delay) that the system exhibited a bistable set of 
solutions. In the rotating state solution, all particles collapse 
to a point and that cluster of particles rotates around a fixed 
center in a circular orbit. The other possible stable solution 
is a ring state, in which all particles distribute themselves 
uniformly along a circle and orbit around its center at unit 
speed. Interestingly, not all particles traverse the ring in 
the same direction; roughly half move clockwise and half 
move anti-clockwise. We will now examine these two states, 
but with random delays given by the truncated Gaussian 
distribution in Eq. ([T3~l >. 

Figure [2] shows the two final particle distributions after 
transients in a simulation with an initial state of N = 150 
randomly placed particles, and where \i x = 2 and c T = 
0.15. In this case, depending upon the random selection of 
delays, either stable solution (ring or rotating) is possible. To 
understand the effects of increasing the standard deviation of 
the random delays, we use a Monte Carlo method. At 100 
different values of a z in the range < (7 T < 0.5, we generate 
random time delays from the distribution in Eq. (TLTt . we then 
simulate the system starting from the same initial condition 
and we determine what state is acquired by the swarm in 
the long-time limit. To determine this, we first measure the 
time-averaged distance of particle j to the center of mass 
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Fig. 3. As (J, is increased, we see a bifurcation from the stable rotating 
state. Panel (a) captures the transition from the rotating state to the ring 
state as the standard deviation of the random delay increases. Panel (b) 
shows the probability of converging to the ring state for a given o~ T of the 
delays. These results were compiled using a Monte Carlo simulation with 
100 random distributions of delays for 100 uniformly-spaced values of <J T 
and for N = 50 particles. See accompanying online movie and Appendix to 
see the agents converge to each stable pattern. 

over the interval (ti, ?2) : 

(8rj) {tl , t2) = -^— f 2 \Srj(t)\dt (17) 

<2 — 'I Jti 

where the size of the interval (fj, ti) is long enough to 
include several periods of oscillation. The ensemble average 
of Eq. ( fT7] i is then 

A value (SrW^) ~ 1 j ' ^fa will indicat^H that the system has 
converged to the ring state, while (8r)r tlitz \ <§; 1/y/a shows 
that the rotating state has been adopted insteacQ. 

Figure |3(a)| demonstrates the effect of increasing <7 T on 
the final state. The blue circles show that for a z small, this 
initial condition converges to the rotating state. However, for 
(7 X > 0.2 the same initial collection of particles will converge 
to the ring state with high probability. In between, there is a 
transition region where both states are commonly observed; 
the state that occurs depends on the random choice of time 
delays. The black dashed lines of |3(a)| show two simulations, 
one which starts near the rotating state (the lower curve), and 
one which starts near the ring state (the upper curve) as <7 T 
is increased. These curves demonstrate the stability of these 
steady states, and the effect of random delays near these 
states. 

Figure |3(b)| shows the conditional probability of ending 
up in the ring state as a function of a T . As expected, for 
this choice of initial conditions, for o T small enough, there 
is zero probability of leaving the rotating state; however, as 
cr T is increased, the probability increases to one. 

The results of these numerical studies strongly suggest 
that even though there is bi-stability between the ring and 
the rotating states, the size of their respective basins of 
attraction is changing dramatically as the standard deviation 
<T T increases. 

'When the delays are uniform, the ring state has a radius of [32]. 
2 This is true for the range of values of o" t considered here. 



VI. Discussion 

In this paper we studied the dynamics of a self-propelling 
swarm with time-delayed inter-agent attraction. In contrast to 
the previously considered case of uniform time delay across 
agents, we considered the situation in which the time delay 
between every pair of agents is drawn randomly from a 
distribution p T . 

Using a mean-field model of the swarm, we showed how 
the two parameter bifurcation plane of coupling strength 
and mean time delay changes with respect to the case in 
which all time delays are equal. The full implications of 
these bifurcation results are the subject of our ongoing work. 
In particular, it is unclear what the stable solutions are. 
Nevertheless, the dramatic changes seen in the two parameter 
bifurcation plane as the standard deviation <7 T increases 
suggest that the basins of attraction of each attractor undergo 
big changes as well. 

Our numerical experiments show that the swarm displays 
bi-stable behavior between the ring and rotating states, at 
the parameters considered. Interestingly, however, our work 
suggests that the basin of attraction of the ring state greatly 
expands as the distribution of time-delays p T widens. Thus, 
in a sense, widening the distribution of time-delays stabilizes 
the stationary state of the swarm center of mass. 

Even though in our model the attractive force among 
agents is linear, we believe this work is useful since it 
represents a first approximation for other, more general forms 
of attractive interaction. Here, we have limited our focus to 
the case where the delays between agents are symmetric and 
constant. However, one important generalization of this sys- 
tem involves incorporating time dependent delays, including 
those which vary as a function of the distance between the 
two agents. This particular refinement of our model is the 
subject of ongoing work and beyond the scope of the current 
paper. 

Finally, although we did not consider repulsion between 
agents, preliminary research leads us to believe that the 
patterns observed in this investigation persist when the 



characteristic repulsion strength between robots is small 
compared to global attraction parameters. For these reasons, 
our results indicate how to exploit time-delayed actuation 
when designing swarm robotic systems with desired tasks 
and functionalities. 

VII. Appendix- Video Description 

The purpose of this research is to investigate the effects 
of randomized communication delay on emerging patterns in 
swarming dynamics. This short video captures the transition 
between two different stable patterns for a swarm as a 
function of the standard distribution of the delays. 

The two coordinate axes in the video show a scatter plot 
of the positions of the particles animated in time. The initial 
positions are identically randomly distributed particles in the 
unit box. The temporal state of the swarm is updated in 
time using a numerical scheme called Heun's Method, and a 
snapshot is captured at every discrete time interval. Here, the 
left coordinate axis uses a standard deviation of the delays 
a x =0.1 while the right axis uses a standard deviation of 
(7 T = 0.3. The mean delay for each simulation is set to be 
jU T = 2, and the number of particles for both is N = 50. 
The vectors at each particle give the velocity associated 
with that particle. The two simulations are run side by side 
to demonstrate the dynamics involved in converging to the 
"rotating" final state on the left, and the "ring" final state 
on the right. The video demonstrates the dynamics over the 
time interval from t = to t = 45, and so includes transients. 
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